

/*
 * Class:        GammaAcceptanceRejectionGen
 * Description:  gamma random variate generators using a method that combines
                 acceptance-rejection with acceptance-complement
 * Environment:  Java
 * Software:     SSJ 
 * Copyright (C) 2001  Pierre L'Ecuyer and Universite de Montreal
 * Organization: DIRO, Universite de Montreal
 * @author       
 * @since

 * SSJ is free software: you can redistribute it and/or modify it under
 * the terms of the GNU General Public License (GPL) as published by the
 * Free Software Foundation, either version 3 of the License, or
 * any later version.

 * SSJ is distributed in the hope that it will be useful,
 * but WITHOUT ANY WARRANTY; without even the implied warranty of
 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
 * GNU General Public License for more details.

 * A copy of the GNU General Public License is available at
   <a href="http://www.gnu.org/licenses">GPL licence site</a>.
 */

package umontreal.iro.lecuyer.randvar;
import umontreal.iro.lecuyer.rng.*;
import umontreal.iro.lecuyer.probdist.*;


/**
 * This class implements <EM>gamma</EM> random variate generators using
 * a method that combines acceptance-rejection 
 * with acceptance-complement.
 * It uses acceptance-rejection for <SPAN CLASS="MATH"><I>&#945;</I> &lt; 1</SPAN> and acceptance-complement
 * for 
 * <SPAN CLASS="MATH"><I>&#945;</I>&nbsp;&gt;=&nbsp;1</SPAN>.
 * For each gamma variate, the first uniform required is taken from the 
 * main stream and all additional uniforms (after the first rejection)
 * are obtained from the auxiliary stream.
 * 
 */
public class GammaAcceptanceRejectionGen extends GammaGen  {
    
   private RandomStream auxStream;
   private static final double q1 =  0.0416666664;
   private static final double q2 =  0.0208333723;
   private static final double q3 =  0.0079849875;
   private static final double q4 =  0.0015746717;
   private static final double q5 = -0.0003349403;
   private static final double q6 =  0.0003340332;
   private static final double q7 =  0.0006053049;
   private static final double q8 = -0.0004701849;
   private static final double q9 =  0.0001710320;
   private static final double a1 =  0.333333333;
   private static final double a2 = -0.249999949;
   private static final double a3 =  0.199999867;
   private static final double a4 = -0.166677482;
   private static final double a5 =  0.142873973;
   private static final double a6 = -0.124385581;
   private static final double a7 =  0.110368310;
   private static final double a8 = -0.112750886;
   private static final double a9 =  0.104089866;
   private static final double e1 =  1.000000000;
   private static final double e2 =  0.499999994;
   private static final double e3 =  0.166666848;
   private static final double e4 =  0.041664508;
   private static final double e5 =  0.008345522;
   private static final double e6 =  0.001353826;
   private static final double e7 =  0.000247453;

   private int gen;
   // UNURAN parameters for the distribution
   private double beta;
   private double gamma;
   // Generator parameters
   // Acceptance rejection combined with acceptance complement
   private static final int gs = 0;
   private static final int gd = 1;
   private double b;
   private double ss;
   private double s;
   private double d;
   private double r;
   private double q0;
   private double c;
   private double si;




   /**
    * Creates a gamma random variate generator with parameters <SPAN CLASS="MATH"><I>&#945;</I> =</SPAN> 
    *  <TT>alpha</TT> and <SPAN CLASS="MATH"><I>&#955;</I> =</SPAN> <TT>lambda</TT>, using main stream <TT>s</TT> and 
    *   auxiliary stream <TT>aux</TT>.
    *  The auxiliary stream is used when a random number of uniforms
    *  is required for a rejection-type generation method.
    * 
    */
   public GammaAcceptanceRejectionGen (RandomStream s, RandomStream aux,
                                       double alpha, double lambda)  {
      super (s, null);
      auxStream = aux;
      setParams (alpha, lambda);
      beta  = 1.0/lambda;
      gamma = 0.0;
      init ();
   }


   /**
    * Creates a gamma random variate generator with parameters <SPAN CLASS="MATH"><I>&#945;</I> =</SPAN> 
    *  <TT>alpha</TT> and <SPAN CLASS="MATH"><I>&#955;</I> =</SPAN> <TT>lambda</TT>, using stream <TT>s</TT>.
    * 
    */
   public GammaAcceptanceRejectionGen (RandomStream s,
                                       double alpha, double lambda)  {
      this (s, s, alpha, lambda);
   }


   /**
    * Creates a new generator object for the gamma 
    *     distribution <TT>dist</TT>, using main stream <TT>s</TT> and 
    *     auxiliary stream <TT>aux</TT>. 
    *     The auxiliary stream is used when a random number of uniforms
    *     is required for a rejection-type generation method.
    * 
    */
   public GammaAcceptanceRejectionGen (RandomStream s, RandomStream aux, 
                                       GammaDist dist)  {
      super (s, dist);
      auxStream = aux;
      if (dist != null)
         setParams (dist.getAlpha(), dist.getLambda());
      beta  = 1.0/dist.getLambda();
      gamma = 0.0;
      init ();
   }


   /**
    * Creates a new generator object for the gamma
    *     distribution <TT>dist</TT> and  stream <TT>s</TT> for both the main and
    *     auxiliary stream.
    * 
    */
   public GammaAcceptanceRejectionGen (RandomStream s, GammaDist dist)  {
      this (s, s, dist);
   }


   /**
    * Returns the auxiliary stream associated with this object.
    * 
    */
   public RandomStream getAuxStream() {
      return auxStream;
   }


   /**
    * Generates a new gamma variate with parameters
    *  <SPAN CLASS="MATH"><I>&#945;</I> =</SPAN>&nbsp;<TT>alpha</TT> and <SPAN CLASS="MATH"><I>&#955;</I> =</SPAN>&nbsp;<TT>lambda</TT>, using
    *   main stream <TT>s</TT> and auxiliary stream <TT>aux</TT>.
    * 
    */
   public static double nextDouble (RandomStream s, RandomStream aux, 
                                    double alpha, double lambda) {
      int gen = gs;
      double s_ = 0, ss = 0, d = 0, q0 = 0, r = 0, c = 0, si = 0, b = 0;

         // Code taken from UNURAN
         if (alpha < 1.0) {
            gen = gs;
            b = 1.0 + 0.36788794412*alpha;       // Step 1
         }
         else {
            gen = gd;
            // Step 1. Preparations
            ss = alpha - 0.5;
            s_ = Math.sqrt (ss);
            d = 5.656854249 - 12.0*s_;

            // Step 4. Set-up for hat case
            r = 1.0 / alpha;
            q0 = ((((((((q9*r + q8)*r + q7)*r + q6)*r + q5)*r + q4)*
                    r + q3)*r + q2)*r + q1)*r;
            if (alpha > 3.686) {
              if (alpha > 13.022) {
                b = 1.77;
                si = 0.75;
                c = 0.1515/s_;
              }
              else {
                b = 1.654 + 0.0076 * ss;
                si = 1.68/s_ + 0.275;
                c = 0.062/s_ + 0.024;
              }
            }
            else {
              b = 0.463 + s_ - 0.178*ss;
              si = 1.235;
              c = 0.195/s_ - 0.079 + 0.016*s_;
            }
         }
         return acceptanceRejection
             (s, aux, alpha, 1.0/lambda, 0, gen, b, s_, ss, d, r, q0, c, si);
   }
 
    
   public double nextDouble() {
      return acceptanceRejection
       (stream, auxStream, alpha, beta, gamma, gen, b, s, ss, d, r, q0, c, si);
   }

   /**
    * Same as <TT>nextDouble (s, s, alpha, lambda)</TT>.
    * 
    */
   public static double nextDouble (RandomStream s, double alpha, 
                                    double lambda) {
      return nextDouble (s, s, alpha, lambda);
   }


   private static double acceptanceRejection
      (RandomStream stream, RandomStream auxStream,
       double alpha, double beta, double gamma, int gen,
       double b, double s, double ss,
       double d, double r, double q0, double c, double si) {
      // Code taken from UNURAN
      double X, p, U, E;
      double q,sign_U,t,v,w,x;
      switch (gen) {
      case gs:
         while (true) {
            p = b*stream.nextDouble();
            stream = auxStream;
            if (p <= 1.0) {                   // Step 2. Case gds <= 1
               X = Math.exp (Math.log (p)/alpha);
               if (Math.log (stream.nextDouble()) <= -X)
                  break;
            }
            else {                           // Step 3. Case gds > 1
               X = -Math.log ((b - p) / alpha);
               if ( Math.log (stream.nextDouble()) <= ((alpha - 1.0)*Math.log (X)))
                  break;
            }
         }
         break;
      case gd:
        do {

            // Step 2. Normal deviate
            t = NormalDist.inverseF01 (stream.nextDouble());
            stream = auxStream;
            x = s + 0.5*t;
            X = x*x;
            if (t >= 0.)
               break;         // Immediate acceptance

            // Step 3. Uniform random number
            U = stream.nextDouble();
            if (d*U <= t*t*t) 
               break;         // Squeeze acceptance

            // Step 5. Calculation of q
            if (x > 0.) {
               // Step 6.
               v = t/(s + s);
               if (Math.abs (v) > 0.25)
                  q = q0 - s*t + 0.25*t*t + (ss + ss)*Math.log (1. + v);
               else
                  q = q0 + 0.5*t*t*((((((((a9*v + a8)*v + a7)*v + a6)*
                                     v + a5)*v + a4)*v + a3)*v + a2)*v + a1)*v;
               // Step 7. Quotient acceptance
               if (Math.log (1. - U) <= q)
                  break;
            }

            // Step 8. Double exponential deviate t
            while (true) {
               do {
                  E = -Math.log (stream.nextDouble());
                  U = stream.nextDouble();
                  U = U + U - 1.;
                  sign_U = (U > 0) ? 1. : -1.;
                  t = b + (E*si)*sign_U;
               } while (t <= -0.71874483771719);   // Step 9. Rejection of t

               // Step 10. New q(t)
               v = t/(s + s);
               if (Math.abs (v) > 0.25)
                  q = q0 - s*t + 0.25*t*t + (ss + ss)*Math.log (1. + v);
               else
                  q = q0 + 0.5*t*t * ((((((((a9*v + a8)*v + a7)*v + a6)*
                                          v + a5)*v + a4)*v + a3)*v + a2)*v + a1)*v;

               // Step 11.
               if (q <= 0.)
                  continue; 

               if (q > 0.5)
                  w = Math.exp (q) - 1.;
               else
                  w = ((((((e7 * q + e6) * q + e5) * q + e4) * q + e3) * q + e2) *
                         q + e1) * q;

               // Step 12. Hat acceptance
               if ( c * U * sign_U <= w*Math.exp (E - 0.5*t*t)) {
                  x = s + 0.5 * t;
                  X = x * x;
                  break;
               }
            }
         } while (false);
         break;
      default: throw new IllegalStateException();
      }

      return gamma + beta*X;
   }

   private void init() {
      // Code taken from UNURAN
      if (alpha < 1.0) {
         gen = gs;
         b = 1.0 + 0.36788794412*alpha;       // Step 1
      }
      else {
         gen = gd;
         // Step 1. Preparations
         ss = alpha - 0.5;
         s = Math.sqrt (ss);
         d = 5.656854249 - 12.0*s;

         // Step 4. Set-up for hat case
         r = 1.0 / alpha;
         q0 = ((((((((q9*r + q8)*r + q7)*r + q6)*r + q5)*r + q4)*
                 r + q3)*r + q2)*r + q1)*r;
         if (alpha > 3.686) {
           if (alpha > 13.022) {
             b = 1.77;
             si = 0.75;
             c = 0.1515/s;
           }
           else {
             b = 1.654 + 0.0076 * ss;
             si = 1.68/s + 0.275;
             c = 0.062/s + 0.024;
           }
         }
         else {
           b = 0.463 + s - 0.178*ss;
           si = 1.235;
           c = 0.195/s - 0.079 + 0.016*s;
         }
      }
   }

}
